source('../env.R')
It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2462 Columns: 10── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities_summary = communities %>% group_by(city_id) %>% summarise(
regional_pool_size = n(),
urban_pool_size = sum(relative_abundance_proxy > 0)
)
ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))
Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data
community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, locomotory_trait_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
left_join(communities_summary) %>%
left_join(city_points[,c('city_id', 'city_nm')]) %>%
rename(city_name='city_nm') %>%
na.omit() %>%
arrange(city_id)
Rows: 341 Columns: 37── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (37): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_mean, fdiv_sd, mass_fdi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`
community_data_metrics
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))
Rows: 304 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
fetch_normalised_traits = function(required_species_list) {
required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
required_traits$trophic_trait_normalised = normalise(required_traits$trophic_trait, min(required_traits$trophic_trait), max(required_traits$trophic_trait))
required_traits$locomotory_trait_normalised = normalise(required_traits$locomotory_trait, min(required_traits$locomotory_trait), max(required_traits$locomotory_trait))
required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), labels = c('Gape Width', 'Trophic Trait', 'Locomotory Trait', 'Mass'))
traits_normalised_long
}
fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))
Read in our phylogeny
phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')
Load resolve ecoregions
resolve = read_resolve()
Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees
to_species_matrix = function(filtered_communities) {
filtered_communities %>%
dplyr::select(city_id, jetz_species_name) %>%
distinct() %>%
mutate(present = TRUE) %>%
pivot_wider(
names_from = jetz_species_name,
values_from = "present",
values_fill = list(present = F)
) %>%
tibble::column_to_rownames(var='city_id')
}
community_nmds = function(filtered_communities) {
species_matrix = to_species_matrix(filtered_communities)
nmds <- metaMDS(species_matrix, k=2, trymax = 30)
nmds_result = data.frame(vegan::scores(nmds)$sites)
nmds_result$city_id = as.double(rownames(nmds_result))
rownames(nmds_result) = NULL
nmds_result
}
https://www.datacamp.com/tutorial/k-means-clustering-r
scree_plot = function(community_nmds_data) {
# Decide how many clusters to look at
n_clusters <- min(10, nrow(community_nmds_data) - 1)
# Initialize total within sum of squares error: wss
wss <- numeric(n_clusters)
set.seed(123)
# Look over 1 to n possible clusters
for (i in 1:n_clusters) {
# Fit the model: km.out
km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
# Save the within cluster sum of squares
wss[i] <- km.out$tot.withinss
}
# Produce a scree plot
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_point(size = 4) +
geom_line() +
geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
xlab('Number of clusters')
scree_plot
}
cluster_cities = function(city_nmds, cities_community_data, centers) {
set.seed(123)
kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
city_nmds$cluster = kmeans_clusters$cluster
cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}
plot_nmds_clusters = function(cluster_cities) {
cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}
get_presence_cell_width = function(city_cluster_data_metrics) {
10 * length(unique(city_cluster_data_metrics$city_id))
}
get_presence_cell_height = function(city_cluster_data_metrics) {
species = species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
dplyr::select(jetz_species_name) %>%
distinct()
10 * nrow(species)
}
city_metric_height = 30
traits_width = 50
phylo_tree_width = 125
title_height = 3
get_image_height = function(city_cluster_data_metrics) {
get_presence_cell_height(city_cluster_data_metrics) + (2 * city_metric_height) + title_height
}
get_image_width = function(city_cluster_data_metrics) {
get_presence_cell_width(city_cluster_data_metrics) + traits_width + phylo_tree_width
}
plot_city_cluster = function(city_cluster_data_metrics, title) {
species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
dplyr::select(jetz_species_name, city_name, relative_abundance_proxy)
tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
gg_tree = ggtree(tree_cropped)
gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) +
geom_tile(aes(fill=relative_abundance_proxy)) +
scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
theme_minimal() + xlab(NULL) + ylab(NULL) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(fill='Urban Proxy Abundance')
species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
gg_cities_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("FDiv") + ylim(0, 1)
gg_title = ggplot() + labs(title = title) + theme_minimal()
gg_presence_height = get_presence_cell_height(city_cluster_data_metrics)
gg_presence_width = get_presence_cell_width(city_cluster_data_metrics)
gg_presence %>% insert_top(gg_cities_mntd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_fd, height = (city_metric_height / gg_presence_height)) %>% insert_left(gg_tree, width = (phylo_tree_width / gg_presence_width)) %>% insert_right(gg_traits, width = (traits_width / gg_presence_width)) %>% insert_top(gg_title, height = (title_height / gg_presence_height))
}
REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')
nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "San Jose" "Los Angeles" "Concord" "Tijuana" "Bakersfield"
[6] "Fresno" "Sacramento" "Mexicali" "Hermosillo" "Las Vegas"
[11] "Phoenix" "Tucson" "Durango" "Portland" "Chihuahua"
[16] "Aguascalientes" "Seattle" "Ciudad Juárez" "San Luis PotosÃ" "Mexico City"
[21] "Saltillo" "Vancouver" "Salt Lake City" "Albuquerque" "Monterrey"
[26] "Nuevo Laredo" "San Antonio" "Denver" "Austin" "Houston"
[31] "Dallas" "Oklahoma City" "Calgary" "New Orleans" "Kansas City"
[36] "Omaha" "St. Louis" "Bradenton" "Tampa" "Minneapolis [Saint Paul]"
[41] "Atlanta" "Orlando" "Louisville" "Chicago" "Indianapolis"
[46] "Milwaukee"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id))
Run 0 stress 0.1005678
Run 1 stress 0.1017503
Run 2 stress 0.1000046
... New best solution
... Procrustes: rmse 0.007231564 max resid 0.03470212
Run 3 stress 0.1000046
... Procrustes: rmse 0.000004044592 max resid 0.00001391239
... Similar to previous best
Run 4 stress 0.1205431
Run 5 stress 0.1012776
Run 6 stress 0.122422
Run 7 stress 0.1005678
Run 8 stress 0.1000046
... Procrustes: rmse 0.000006352527 max resid 0.00001951675
... Similar to previous best
Run 9 stress 0.1217145
Run 10 stress 0.1000046
... New best solution
... Procrustes: rmse 0.000003781293 max resid 0.00001247782
... Similar to previous best
Run 11 stress 0.1217437
Run 12 stress 0.1217437
Run 13 stress 0.1000046
... Procrustes: rmse 0.000004334204 max resid 0.00001530842
... Similar to previous best
Run 14 stress 0.1000046
... Procrustes: rmse 0.000004734417 max resid 0.00001660964
... Similar to previous best
Run 15 stress 0.1012776
Run 16 stress 0.1235457
Run 17 stress 0.1000046
... Procrustes: rmse 0.000004211501 max resid 0.00001320069
... Similar to previous best
Run 18 stress 0.1005678
Run 19 stress 0.1005678
Run 20 stress 0.1017503
*** Best solution repeated 4 times
nearctic_cities_nmds
scree_plot(nearctic_cities_nmds)
nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 4)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(nearctic_cities)
nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = nearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))
Saving 7.29 x 4.51 in image
nearactic_cluster1 = nearctic_cities %>% filter(cluster == 1)
plot_city_cluster(nearactic_cluster1, 'Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'), width = get_image_width(nearactic_cluster1), height = get_image_height(nearactic_cluster1), units = "mm")
nearactic_cluster2 = nearctic_cities %>% filter(cluster == 2)
plot_city_cluster(nearactic_cluster2, 'Neartic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'), width = get_image_width(nearactic_cluster2), height = get_image_height(nearactic_cluster2), units = "mm")
nearactic_cluster3 = nearctic_cities %>% filter(cluster == 3)
plot_city_cluster(nearactic_cluster3, 'Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'), width = get_image_width(nearactic_cluster3), height = get_image_height(nearactic_cluster3), units = "mm")
nearactic_cluster4 = nearctic_cities %>% filter(cluster == 4)
plot_city_cluster(nearactic_cluster4, 'Neartic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster5.jpg'), width = get_image_width(nearactic_cluster4), height = get_image_height(nearactic_cluster4), units = "mm")
neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Culiacán" "Guadalajara" "Morelia" "Acapulco" "Querétaro"
[6] "Cuernavaca" "Puebla" "Oaxaca" "Xalapa" "Veracruz"
[11] "Tuxtla Gutiérrez" "Villahermosa" "Guatemala City" "San Salvador" "San Pedro Sula"
[16] "Mérida" "Tegucigalpa" "Managua" "San José" "Cancún"
[21] "Guayaquil" "Chiclayo" "Panama City" "Trujillo" "Quito"
[26] "Havana" "Cali" "Lima" "Pereira" "Miami"
[31] "MedellÃn" "Ibagué" "Cartagena" "Kingston" "Bogota"
[36] "Barranquilla" "Bucaramanga" "Cúcuta" "Maracaibo" "Arequipa"
[41] "Barquisimeto" "Santo Domingo" "Maracay" "El Alto [La Paz]" "Caracas"
[46] "Cochabamba" "Viña del Mar [ValparaÃso]" "RÃo Piedras [San Juan]" "Barcelona" "Concepción"
[51] "Santiago" "Mendoza" "Salta" "Cordoba" "Asuncion"
[56] "Buenos Aires" "La Plata" "Ciudad del Este" "Montevideo" "Mar del Plata"
[61] "Porto Alegre" "São Paulo" "Santos" "Sao Jose dos Campos"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id))
Run 0 stress 0.134619
Run 1 stress 0.134433
... New best solution
... Procrustes: rmse 0.006824014 max resid 0.04572879
Run 2 stress 0.1407144
Run 3 stress 0.1348237
... Procrustes: rmse 0.006646789 max resid 0.04604869
Run 4 stress 0.1348046
... Procrustes: rmse 0.006551306 max resid 0.04608394
Run 5 stress 0.134433
... Procrustes: rmse 0.00003293268 max resid 0.0001049247
... Similar to previous best
Run 6 stress 0.134433
... Procrustes: rmse 0.00005189183 max resid 0.0001643412
... Similar to previous best
Run 7 stress 0.1366488
Run 8 stress 0.1407144
Run 9 stress 0.1348046
... Procrustes: rmse 0.006551801 max resid 0.04608656
Run 10 stress 0.134619
... Procrustes: rmse 0.006841563 max resid 0.04571578
Run 11 stress 0.1348237
... Procrustes: rmse 0.006648345 max resid 0.04605669
Run 12 stress 0.1344331
... Procrustes: rmse 0.00006317981 max resid 0.0001865383
... Similar to previous best
Run 13 stress 0.1405637
Run 14 stress 0.1346406
... Procrustes: rmse 0.007270798 max resid 0.04575267
Run 15 stress 0.1405831
Run 16 stress 0.1406945
Run 17 stress 0.1346406
... Procrustes: rmse 0.007285202 max resid 0.04573278
Run 18 stress 0.1344509
... Procrustes: rmse 0.001186371 max resid 0.006935107
... Similar to previous best
Run 19 stress 0.1405634
Run 20 stress 0.1344331
... Procrustes: rmse 0.00007290489 max resid 0.0002305292
... Similar to previous best
*** Best solution repeated 5 times
neotropic_cities_nmds
scree_plot(neotropic_cities_nmds)
neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(neotropic_cities)
neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
ggplot() +
geom_sf(data = neotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cluster1 = neotropic_cities %>% filter(cluster == 1)
plot_city_cluster(neotropic_cluster1, 'Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'), width = get_image_width(neotropic_cluster1), height = get_image_height(neotropic_cluster1), units = "mm")
neotropic_cluster2 = neotropic_cities %>% filter(cluster == 2)
plot_city_cluster(neotropic_cluster2, 'Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'), width = get_image_width(neotropic_cluster2), height = get_image_height(neotropic_cluster2), units = "mm")
neotropic_cluster3 = neotropic_cities %>% filter(cluster == 3)
plot_city_cluster(neotropic_cluster3, 'Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'), width = get_image_width(neotropic_cluster3), height = get_image_height(neotropic_cluster3), units = "mm")
neotropic_cluster4 = neotropic_cities %>% filter(cluster == 4)
plot_city_cluster(neotropic_cluster4, 'Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'), width = get_image_width(neotropic_cluster4), height = get_image_height(neotropic_cluster4), units = "mm")
neotropic_cluster5 = neotropic_cities %>% filter(cluster == 5)
plot_city_cluster(neotropic_cluster5, 'Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'), width = get_image_width(neotropic_cluster5), height = get_image_height(neotropic_cluster5), units = "mm")
palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Lisbon" "Porto" "Marrakesh" "Seville" "Dublin" "Málaga"
[7] "Madrid" "Glasgow" "Bilbao" "Liverpool" "Bristol" "Manchester"
[13] "Birmingham" "Leeds" "Newcastle upon Tyne" "Sheffield" "Nottingham" "Valencia"
[19] "London" "Toulouse" "Paris" "Barcelona" "Rotterdam [The Hague]" "Brussels"
[25] "Amsterdam" "Lyon" "Marseille" "Dusseldorf" "Nice" "Frankfurt am Main"
[31] "Zurich" "Oslo" "Stuttgart" "Hamburg" "Genoa" "Nuremberg"
[37] "Copenhagen" "Munich" "Berlin" "Dresden" "Rome" "Prague"
[43] "Stockholm" "Poznan" "Vienna" "Wroclaw" "Zagreb" "Gdansk"
[49] "Budapest" "Krakow" "Warsaw" "Helsinki" "Riga" "Belgrade"
[55] "Lviv" "Sofia" "Thessaloniki" "Saint Petersburg" "Minsk" "Athens"
[61] "Kyiv" "Istanbul" "Odesa" "Samsun" "Luxor" "Tel Aviv"
[67] "Jerusalem" "Tbilisi" "Yerevan" "Kuwait City" "Doha" "Abu Dhabi"
[73] "Dubai" "Bishkek"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id))
Run 0 stress 0.04961857
Run 1 stress 0.06522066
Run 2 stress 0.06770805
Run 3 stress 0.06505779
Run 4 stress 0.08915955
Run 5 stress 0.05395855
Run 6 stress 0.04961847
... New best solution
... Procrustes: rmse 0.0001353415 max resid 0.0006424586
... Similar to previous best
Run 7 stress 0.07723898
Run 8 stress 0.06982598
Run 9 stress 0.05002297
... Procrustes: rmse 0.008523462 max resid 0.02080155
Run 10 stress 0.05001003
... Procrustes: rmse 0.06742951 max resid 0.2228023
Run 11 stress 0.05354024
Run 12 stress 0.07364969
Run 13 stress 0.04967779
... Procrustes: rmse 0.02454439 max resid 0.1155996
Run 14 stress 0.06518621
Run 15 stress 0.0547908
Run 16 stress 0.0675422
Run 17 stress 0.06304729
Run 18 stress 0.04968731
... Procrustes: rmse 0.02483805 max resid 0.1165918
Run 19 stress 0.06829409
Run 20 stress 0.08070857
*** Best solution repeated 1 times
palearctic_cities_nmds
scree_plot(palearctic_cities_nmds)
palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 7)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(palearctic_cities)
palearctic_biomes = resolve[resolve$REALM == 'Palearctic',c('REALM')]
ggplot() +
geom_sf(data = palearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cluster1 = palearctic_cities %>% filter(cluster == 1)
plot_city_cluster(palearctic_cluster1, 'Palearctic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'), width = get_image_width(palearctic_cluster1), height = get_image_height(palearctic_cluster1), units = "mm")
palearctic_cluster2 = palearctic_cities %>% filter(cluster == 2)
plot_city_cluster(palearctic_cluster2, 'Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'), width = get_image_width(palearctic_cluster2), height = get_image_height(palearctic_cluster2), units = "mm")
palearctic_cluster3 = palearctic_cities %>% filter(cluster == 3)
plot_city_cluster(palearctic_cluster3, 'Palearctic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'), width = get_image_width(palearctic_cluster3), height = get_image_height(palearctic_cluster3), units = "mm")
palearctic_cluster4 = palearctic_cities %>% filter(cluster == 4)
plot_city_cluster(palearctic_cluster4, 'Palearctic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'), width = get_image_width(palearctic_cluster4), height = get_image_height(palearctic_cluster4), units = "mm")
palearctic_cluster5 = palearctic_cities %>% filter(cluster == 5)
plot_city_cluster(palearctic_cluster5, 'Palearctic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")
palearctic_cluster6 = palearctic_cities %>% filter(cluster == 6)
plot_city_cluster(palearctic_cluster6, 'Palearctic cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'), width = get_image_width(palearctic_cluster6), height = get_image_height(palearctic_cluster6), units = "mm")
palearctic_cluster7 = palearctic_cities %>% filter(cluster == 7)
plot_city_cluster(palearctic_cluster7, 'Palearctic cluster 7')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster7.jpg'), width = get_image_width(palearctic_cluster7), height = get_image_height(palearctic_cluster7), units = "mm")
afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Cape Town" "Johannesburg" "Pretoria" "Kigali" "Kampala" "Arusha" "Nairobi" "Addis Ababa" "Antananarivo"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id))
Run 0 stress 0.00009014786
Run 1 stress 0.0000894423
... New best solution
... Procrustes: rmse 0.0002302813 max resid 0.0003204519
... Similar to previous best
Run 2 stress 0.0002552434
... Procrustes: rmse 0.001471005 max resid 0.001947622
... Similar to previous best
Run 3 stress 0.00009602873
... Procrustes: rmse 0.0003134658 max resid 0.0004532654
... Similar to previous best
Run 4 stress 0.0004206102
... Procrustes: rmse 0.002382671 max resid 0.003350133
... Similar to previous best
Run 5 stress 0.001267824
Run 6 stress 0.00009984846
... Procrustes: rmse 0.000573798 max resid 0.0007724523
... Similar to previous best
Run 7 stress 0.00009404138
... Procrustes: rmse 0.0002198999 max resid 0.0002986639
... Similar to previous best
Run 8 stress 0.00006937968
... New best solution
... Procrustes: rmse 0.0001548829 max resid 0.0003491895
... Similar to previous best
Run 9 stress 0.00009515083
... Procrustes: rmse 0.0002518215 max resid 0.0003915592
... Similar to previous best
Run 10 stress 0.0002810887
... Procrustes: rmse 0.001624066 max resid 0.002239305
... Similar to previous best
Run 11 stress 0.3209238
Run 12 stress 0.0001607375
... Procrustes: rmse 0.0009389978 max resid 0.001284128
... Similar to previous best
Run 13 stress 0.00009851387
... Procrustes: rmse 0.0002528255 max resid 0.0003960776
... Similar to previous best
Run 14 stress 0.00009299962
... Procrustes: rmse 0.0002384291 max resid 0.0003789721
... Similar to previous best
Run 15 stress 0.0001360544
... Procrustes: rmse 0.0007908028 max resid 0.001108713
... Similar to previous best
Run 16 stress 0.00009817612
... Procrustes: rmse 0.0002561838 max resid 0.0003936572
... Similar to previous best
Run 17 stress 0.002120919
Run 18 stress 0.00009438058
... Procrustes: rmse 0.000245643 max resid 0.0003771271
... Similar to previous best
Run 19 stress 0.00009470345
... Procrustes: rmse 0.0002746505 max resid 0.0004407572
... Similar to previous best
Run 20 stress 0.00008594013
... Procrustes: rmse 0.0002324527 max resid 0.0003605282
... Similar to previous best
*** Best solution repeated 11 times
Warning: stress is (nearly) zero: you may have insufficient data
afrotropic_cities_nmds
scree_plot(afrotropic_cities_nmds)
afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 2)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(afrotropic_cities)
afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
ggplot() +
geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cluster1 = afrotropic_cities %>% filter(cluster == 1)
plot_city_cluster(afrotropic_cluster1, 'Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'), width = get_image_width(afrotropic_cluster1), height = get_image_height(afrotropic_cluster1), units = "mm")
afrotropic_cluster2 = afrotropic_cities %>% filter(cluster == 2)
plot_city_cluster(afrotropic_cluster2, 'Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'), width = get_image_width(afrotropic_cluster2), height = get_image_height(afrotropic_cluster2), units = "mm")
indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Srinagar" "Jamnagar" "Jammu" "Rajkot" "Bikaner" "Jodhpur" "Jalandhar"
[8] "Ahmedabad" "Bhavnagar" "Ludhiana" "Anand" "Udaipur" "Surat" "Vadodara"
[15] "Ajmer" "Chandigarh" "Vasai-Virar" "Mumbai" "Jaipur" "Delhi [New Delhi]" "Nashik"
[22] "Dehradun" "Kota" "Pune" "Haridwar" "Dhule" "Ujjain" "Indore"
[29] "Ahmadnagar" "Kolhapur" "Jalgaon" "Agra" "Aurangabad" "Sangli" "Belagavi"
[36] "Gwalior" "Budaun" "Bareilly" "Dharwad" "Bhopal" "Bhind" "Mangaluru"
[43] "Solapur" "Vijayapura" "Akola" "Latur" "Kannur" "Davanagere" "Thalassery"
[50] "Amravati" "Kalaburagi" "Kozhikode" "Guruvayur" "Malappuram" "Lucknow" "Thrissur"
[57] "Mysuru" "Kochi" "Alappuzha" "Nagpur" "Kollam" "Jabalpur" "Ettumanoor"
[64] "Hyderabad" "Coimbatore" "Bengaluru" "Thiruvananthapuram" "Tiruppur" "Faizabad" "Erode"
[71] "Prayagraj" "Pratapgarh" "Salem" "Dindigul" "Madurai" "Tiruchirappalli" "Durg"
[78] "Vellore" "Tirupati" "Raipur" "Bilaspur" "Vijayawada" "Puducherry" "Chennai"
[85] "Kathmandu" "Colombo" "Rajamahendravaram" "Patna" "Kandy" "Bihar Sharif" "Visakhapatnam"
[92] "Ranchi" "Brahmapur" "Jamshedpur" "Darjeeling" "Siliguri" "Cuttack" "Bhubaneshwar"
[99] "Jalpaiguri" "Berhampore" "Kolkata" "Krishnanagar" "Guwahati [Dispur]" "Agartala" "Silchar"
[106] "Dimapur" "Bangkok" "George Town" "Kuala Lumpur" "Phnom Penh" "Singapore" "Hong Kong"
[113] "Sha Tin" "Hsinchu" "Taichung" "New Taipei [Taipei]" "Tainan" "Denpasar" "Kaohsiung"
[120] "Kota Kinabalu"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id))
Run 0 stress 0.1190668
Run 1 stress 0.1224401
Run 2 stress 0.1394961
Run 3 stress 0.1161598
... New best solution
... Procrustes: rmse 0.02798297 max resid 0.2346031
Run 4 stress 0.1540252
Run 5 stress 0.1175538
Run 6 stress 0.1167657
Run 7 stress 0.1384412
Run 8 stress 0.1241275
Run 9 stress 0.1163498
... Procrustes: rmse 0.025336 max resid 0.2481842
Run 10 stress 0.1153501
... New best solution
... Procrustes: rmse 0.008287738 max resid 0.08217213
Run 11 stress 0.1199235
Run 12 stress 0.1172201
Run 13 stress 0.151756
Run 14 stress 0.1246307
Run 15 stress 0.134023
Run 16 stress 0.1580175
Run 17 stress 0.1229198
Run 18 stress 0.1502633
Run 19 stress 0.1189366
Run 20 stress 0.117072
Run 21 stress 0.1458634
Run 22 stress 0.137476
Run 23 stress 0.1186188
Run 24 stress 0.1637221
Run 25 stress 0.1191928
Run 26 stress 0.1164559
Run 27 stress 0.142366
Run 28 stress 0.1379337
Run 29 stress 0.1191853
Run 30 stress 0.1266232
*** Best solution was not repeated -- monoMDS stopping criteria:
29: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
indomalayan_cities_nmds
scree_plot(indomalayan_cities_nmds)
indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(indomalayan_cities)
indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
ggplot() +
geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) +
geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cluster1 = indomalayan_cities %>% filter(cluster == 1)
plot_city_cluster(afrotropic_cluster2, 'Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'), width = get_image_width(indomalayan_cluster1), height = get_image_height(indomalayan_cluster1), units = "mm")
indomalayan_cluster2 = indomalayan_cities %>% filter(cluster == 2)
plot_city_cluster(afrotropic_cluster2, 'Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'), width = get_image_width(indomalayan_cluster2), height = get_image_height(indomalayan_cluster2), units = "mm")
indomalayan_cluster3 = indomalayan_cities %>% filter(cluster == 3)
plot_city_cluster(afrotropic_cluster2, 'Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'), width = get_image_width(indomalayan_cluster3), height = get_image_height(indomalayan_cluster3), units = "mm")
indomalayan_cluster4 = indomalayan_cities %>% filter(cluster == 4)
plot_city_cluster(afrotropic_cluster2, 'Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'), width = get_image_width(indomalayan_cluster4), height = get_image_height(indomalayan_cluster4), units = "mm")
indomalayan_cluster5 = indomalayan_cities %>% filter(cluster == 5)
plot_city_cluster(afrotropic_cluster2, 'Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'), width = get_image_width(indomalayan_cluster5), height = get_image_height(indomalayan_cluster5), units = "mm")